Elongation dynamics of amyloid fibrils: a rugged energy landscape picture 
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Protein amyloid fibrils are a form of linear protein aggregates that are implicated in many neu- 
rodegenerative diseases. Here, we study the dynamics of amyloid fibril elongation by performing 
Langevin dynamic simulations on a coarse-grained model of peptides. Our simulation results sug- 
gest that the elongation process is dominated by a series of local minimum due to frustration in 
monomer- fibril interactions. This rugged energy landscape picture indicates that the amount of re- 
cycling of monomers at the fibrils' ends before being fibrilized is substantially reduced in comparison 
to the conventional two-step elongation model. This picture, along with other predictions discussed, 
can be tested with current experimental techniques. 

PACS numbers: 87.14.ef, 82.35.Pq, 83.10.Mj, 46.25.Cc 



I. INTRODUCTION 

Amyloids are insoluble fibrous protein aggregations 
stabilized by a network of hydrogen bonds and hydropho- 
bic interactions [l|, [H, 0, • They are intimately related 
to many neurodegenerative diseases such as Alzheimer's 
Disease, Parkinson's Disease and prion diseases [5]. Bet- 
ter characterization of the various properties of amyloid 
fibrils is therefore of high importance for the understand- 
ing of the associated pathogenesis. In this paper, we in- 
vestigate the dynamics of elongation in fibril growth. 

From the kinetic theory point of view, the elonga- 
tion process is traditionally viewed as a diffusion-limited 
reaction coupled with high free-energy barrier crossing 
m|. A more refined picture has also been proposed in 
B S lid, lUl, Il2l where the elongation process is 
treated as a two-step dock- and- convert process (c.f. Fig. 
[TJa)). Namely, the elongation process follows the kinetics 
scheme: 

[m] + [h]^[{m®fk)]-[fk+i\ (1) 
a_ 

where [m] denotes the monomer concentration, [fk] the 
fibril concentration consisting of k monomers, [(m fk)] 
denotes the intermediate state before a monomer can 
be converted into fibril, and a+, a_ and h are the 
corresponding rates. Experimentally, elongation rates 
for amyloid fibrils formed from various peptides have 
been measured (e.g., [TO, In particular, 
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FIG. 1: (a) The dock- and- convert free energy landscape pic- 
ture of the elongation process proposed in 0, 0, [lol , |Tl[. 
5*0 denotes the initial state, i.e., a free monomer, and S'5 de- 
notes the final state with monomer being part of the fibril, 
(b) The free energy landscape picture advocated in this work. 
The schematic underneath the landscape picture depicts one 
particular conformation in the 6*2^3 state trapped due to mis- 
alignment. Note that there are three D-bonds formed with 
the fibril and the cloud encloses the dangling end with two 
unbound A beads. 



(a) m 




of Alzheimer's disease [5|, have an estimated conversion- 
limiting elongation rate of ~0.3 /im/min [14]. Given that 
each beta strand is about 0.5 nm in width, this elonga- 
tion rate translates to be in the order of one monomer 
fibrilized per second, i.e., 6+ ~ 1 s~^. Intuitively, the in- 
termediate state {m® fk) corresponds to the initial state 
when the monomer first becomes bound to the fibril's end 
through hydrophobic interactions or /and a small number 
of hydrogen bonds. This suggests that at the initial bind- 
ing, the energy gain, A/7, should amount to only a few 
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ksT fid\. If we estimate a_ by treating the undocking 
process as a diffusion-hmited dissociation , then (see, e.g., 
ch. 8 in [M]): 



3De 



AU 



(2) 



where r is the distance between the two react ants, which 
is of the order of 1 nm. Taking the binding energy, AU to 
be —SksT^ say, and setting D to be 10~^^ m^s~^, which 
is a typical diffusion constant for small proteins [l^, a_ 
can be estimated to be about 10'' s~^. Since 6+ ~ 1 s~^ 
and a_ ~ 10^ s~^, the two-step model depicted in Eq. 
([1]) suggests that almost 10^ monomers would have been 
interacted with the fibril's end before one of them is con- 
verted [13 under the experimental condition investigated 
in 

In this work, we demonstrate, with the help of coarse- 
grained molecular dynamics simulations, that the dy- 
namics in the conversion step is dominated by a series 
of energy traps manifested by the frustrated monomer- 
fibril interactions, which would include i) the misalign- 
ment of the monomer with respect to the beta strands 
at the fibril's end; ii) frustrated hydrophobic interactions 
among the side chains; and iii) the competition to bind 
between multiple monomers at the same end of a grow- 
ing fibril. In this picture, there is not a rate limiting step 
for the elongation process, but rather, each energy trap 
contributes to the final elongation rate observed. This 
scenario is akin to the random energy model investigated 
in glassy systems (e.g., see [Hi]). This rugged energy 
landscape picture predicts that monomers would spend 
a substantial amount of time at the fibril's end before 
conversion. As a result, the amount of recycling of the 
monomers at the fibrils' ends before one of them becomes 
fibrilized would be many orders of magnitude small than 
the value indicated by the two-step model depicted in Eq. 
[TJ This dynamical picture is testable by, e.g., perform- 
ing fibril elongation experiments with a small portion of 
monomers radioactively tagged [10|, 0]. 

We will now present the details on the simulation 
method is presented in Section [TTl In Section Hill we 
explain how the simulation results are analyzed. In Sec- 
tion [iVl we discuss how our findings relate to the rugged 
energy picture introduced and present some implications 
of our model. We then end with a conclusion. 



II. SIMULATION DETAILS 

Employing coarse-grained peptide models for the study 
of amyloids are abundant in the literature (e.g., [H, [H, 
HqI, El, S m, m, [13), here we aim to use the simplest 
model to study amyloid fibril elongation. Specifically, 
we keep only two types of inter-peptide interactions: di- 
rectional interactions provided by hydrogen bonds, and 
undirectional interactions given by hydrophobic interac- 
tions. The directional interactions dictate that fibrils can 
only grow linearly, and the undirectional interactions are 



FIG. 2: Schematic pictures depicting the model adopted in 
this work, (a) Each amino acid is simplistically represented 
by two beads, the gray beads represent the peptide back- 
bone of an amino acid and the coloured beads the side chains 
(red for hydrophobic and green for hydrophilic) . We stress 
that the representation is only meant to be qualitative, (b) 
The five- amino- acid peptide employed in this work with alter- 
nating hydrophilic-hydrophobic side chains. The alternating 

P' tern has been shown to promote amyloid fibril formation 
, [13, M, [H, [13, m. (c) A segment of fibril consisting of 
four peptides in two layers of cross-beta sheets, (d) A car- 
toon depicting the four beta strands corresponding to (c). 
The green (red) face of the panel depicts the hydrophilic (hy- 
drophobic) side. 



(a) Si| 




(b) 




(d) 




indispensable in keeping the fibrils thermodynamically 
stable [26] . The directionality of the hydrogen bonds sug- 
gests that each amino acid has to be represented by at 
least two set of coordinates, one for its position and one 
for the direction of the side-chain. We therefore min- 
imally employ two beads to represent each amino-acid 
(c.f. Fig.lSfa)). We stress that we do not attempt to de- 
vise a quantitatively correct representation of a peptide, 
but rather to use a toy model to study the elongation 
process. 

We will now describe the details of the model. For 
simplicity, there are only two types of interactions: 



F(r, hz^ cr, d) 
W{r, cr, d) 







\r - a\ < d 
otherwise 

, r < d 
otherwise 



(3) 



(4) 



Namely, F is a harmonic potential with a sharp cutoff at 
d (in units of nm) and W is the Lennard- Jones poten- 
tial again with a cutoff at a (in units nm) denotes the 
minimum in the potential well and controls the depth 
of the potential well (in units of ksT). Note the dis- 
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FIG. 3: The bonded interactions for the model, (a) Besides 
the white bonds between the beads, the orange bonds are 
employed to fix the angle between the side-chains and the 
peptide backbone. Note that the side-chain at the end of the 
peptide, B5, is bonded to A4 in order to maintain the angle 
B5—A5 — A4. (b) The six bonding interactions in the D-bond 
between beads Ak and Ah. If one of the A beads involved is 
at the end of the peptide, e.g.. Ah = A5, then the cross- 
peptide A- A bonds will be between are Ak — A5, Afc+i — A5 
and Ak — A4 instead. 




continuities in the slopes of the potentials at the cutoffs. 
We believe that these discontinuities are unimportant in 
our Langevin simulations due to the greater magnitude 
of perturbation from thermal fluctuation . 

The whole system is modeled by pairwise interactions 
consisting of a linear sum of a set of potentials V and 
which can be categorized into two classes: bonded 
interactions and non-bonded interactions. 



TABLE I: (a) The parameters employed in the bonded inter- 
actions. The potential is of the form V with d = 00 and the 
values in the entries correspond to {k, a), (b) The parameters 
for [/-type non-bonded interactions. The potential is of the 
form W and the values in the entries correspond to {k, a, d). 
(c) The parameters for D-type non-bonded interactions. The 
potential is of the form V and the values in the entries corre- 
spond to {K,a,d). Note that these interaction potentials are 
only switched on when all six distances are within the cutoffs. 



(a) Bonded interactions 
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(b) ^7-type non-bonded interactions between 
every pair of beads not bonded 
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(5,1/2,1/2) 


^i=odd 


(5,1/2,1/2) 


(5,1/2,1/2) 


(3/2,1/2,1) 



(c) D-type non-bonded interactions between 
beads and^j 





A. 








(4,1/2,1/10) 


(4,2-1^2^3/20) 


(4,2-1^2^3/20) 




(4,2-1^2^3/20) 


(4,1/2,1/10) 




^1+1 


(4,2-1^2^3/20) 







Non-bonded interactions 



A. Bonded interactions 

Bonded interactions refer to the bonds within a pep- 
tide in order to provide it with the structural constraints 
that mimic a peptide (c.f. Fig. [3]). The parameters of 
the interactions are shown in Table 1(a). For example, 
the total force acting on the bead A\ due to bonded in- 
teractions is: 

Va, \y{TA,A^ , 40, 1/2, 00) + V{ta,b, , 40, 1/2, 00) 
+ F(rA,A3, 10,1,00)1 , (5) 



where = Vr^ and Tol(5 = Itq, — r/3| with y^, referring 
to the position of the a bead. The first term fixes the 
inter- amino- acid distance to be about 0.5 nm and the 
second term dictates that the distance between the pep- 
tide backbone and the side chain to be about 0.5 nm. 
The second term above gives longitudinal rigidity to the 
peptide. 



Besides the bonded interactions, there are also non- 
bonded interactions between the beads. The first type is 
undirectional and we call them [/-type interactions. 



U-type interactions 



These interactions include steric constraints or at- 
tractive interactions (only between pairs of hydrophobic 
beads, i.e., 5i=odd) between every pair of beads, except 
for pairs already under bonded interactions. These ef- 
fects are manifested by the the Lennard- Jones potential 
with different cutoffs, d: d = a for the case of pure re- 
pulsion, and d > a for the case of long range attraction 
with short range repulsion. The parameters employed in 
the simulations are shown in Table lb. For instance, the 
interaction potential between two red beads, Bi and Bj^ 
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2. D-type interactions 

We use the term I)-bonds to refer to the directional 
interactions between peptides. The I)-bonds are meant 
to mimic the cross-beta sheet hydrogen bonds. We will 
model the directional elements by a sum of six harmonic 
potentials V (c.f. Fig.[3l^b)). This way of modeling direc- 
tionality in hydrogen bonding is akin to the works em- 
ploying discrete molecular dynamics simulations to study 
amyloid formation in the literature [19]. Note that the 
potentials are only switched on when all six pairwise dis- 
tances concerned are within their respective cutoffs. In 
other words, the total interaction potential for a D-type 
bond between Ai and Aj is: 



VivA^A, , 5, 1/2, 1/10) + V{rA^B, , 5, 2-l/^ 3/20) 
V{rA^A,+, , 5, 2-1/2, 3/20) + V{rB^A, , 5, 2''/^ 3/20) 
V{rB^B, , 5, 1/2, 1/10) + V{rA^^,A, , 5, 2-^/^, 3/20) 



TABLE II: A comparison between the parameters employed 
in the simulations measure (middle column) and the cor- 
responding values measured experimentally (right column). 
Note that the hydrophobic strength and the hydrogen bond 
strength in the middle column depicts the enthalpies defined 
in the simulations while the values in the right columns are 
measured in free energies. We note that the qualitative na- 
ture of our conclusion stays the same when the hydrophobic 
strength (hydrogen bond strength) are varied around the pre- 
sented values by twenty (ten) percents above and below. 



Properties 


Simul. 


Exp'l 


Friction coefficient per bead, 7 (/c^Tps/nm^) 


1000 


- 1000 [27] 


Inter- amino- acid distance (nm) 


0.5 


0.35 [28] 


Hydrophobic interaction strength (ksT) 


1.5 


1-4 [16] 


Hydrogen bond strength {ksT) 


4 


- 2.3 [16] 


Time increment, At (fs) 


5.6 





where 9 = 1 when all of the six distances are within their 
respective cutoffs, and 6> = otherwise. 



C. Simulation procedure 

We performed Langevin dynamics simulations on our 
system. Namely, the position of the a bead, r^, follows 
the updating rule 0: 



r«(t+At) = r«(t) 



-V«/7totai({r})At+W?^At7f« , 
7 V ^ 

(6) 

where t/q, represents Gaussian noise with zero mean and 
variance one in 3D, 7 is the friction coefficient for each 
bead, and /7totai is the sum of all pairwise interactions in 
the model. The relevant parameters are shown in Table 
II. 

Simulations are done with one fibril segment and one 
monomer in a cubic box 6 nm on a side (therefore, the 
monomer concentration [m] and fibril concentration [/] 
is 7.7 mM). A fibril segment consisting of ten 5-amino- 
acid peptides are placed at the center of the box. The 
fibril is constructed by hand and consists of two-layer 
of cross-beta sheet structure as depicted in Fig. [H The 
fibril is held fixed, i.e., the peptides within it are com- 
pletely frozen throughout the simulation. At time zero, 
a monomeric peptide is placed at the corner of the box 
and the simulation is stopped when the free monomer has 
all of the five I)-bonds formed with the fibril. Note that 
there are four possible locations for the added monomer 
to bind to as the fibril has two ends and there are two 
cross-beta sheets. 

Throughout the run, we record the time when a change 
in the number of I)-bonds between the monomeric pep- 
tide and the fibril happens. This allows us to construct a 



FIG. 4: A snapshot of one simulation run in progress. 
The fibril is put in the middle of the simulation box and 
the monomer is diffusing towards it. The fibrillar axis is 
along the 2;-axis, the cross-beta sheets are along the x-axis, 
and the coordinates of the A^ bead for the peptides are 
(-0.25, 0.6, k/2 - 1) and (0.25, 0.6, k/2 - 1.25), /c = 0, . . . 4. 
The corners of the simulation box are at x,y,z = =b3. 




time series describing the temporal evolution of the elon- 
gation process. We will now describe how the time series 
is analyzed. 
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III. DATA ANALYSIS 

To comprehend our simulation results, we adopt a 
coarse-grained picture of the dynamics of fibril elonga- 
tion. Specifically, we partition the phase space of a 
monomer in the process of fibrilization into a number 
of discrete states, and aim to approximate the dynamical 
picture by a series of jumps between neighboring states 
by Markovian processes. The desire to have a Marko- 
vian representation is an attempt to view the dynamics 
through a familiar mechanism. To find a sensible def- 
inition for the set of discrete states, we firstly simplify 
the dynamical picture by recoding the time whenever a 
i^-bond is formed or destroyed. This gives an array con- 
sisting of the times and the numbers of I)-bonds between 
the monomer and the fibril as shown in Fig. [Sja). By 
inspection of the dataset, it is apparent that the time 
series consists of segments of long periods within which 
there are a lot of rapid back and forth transitions be- 
tween having k and k -\- 1 I)-bonds. This is a clear sign 
of temporal correlation and since our desire is to approx- 
imate the process with a memoryless kinetic mechanism, 
we will partition the configuration space of our system 
into five discrete states designated by: So, ^1^2, ^2^3, 
Ss^4 and ^5, where refers to having no A-bonds be- 
tween the monomer and the fibril, Sk^k-\-i refers to the 
state where the number of I)-bonds fiickers between k 
and /c + 1, and ^5 refers to the fully aligned state for the 
monomer. With these newly defined states, a new time 
series recording the transitions between them can be con- 
structed (c.f. Fig. [5] and Fig. [6]). We now assume that 
all the transition events are drawn from Independent and 
Exponential Distributions. Given the property that the 
minimum of two exponential random variables is again 
exponentially distributed with rate equal to the sum of 
the two original rates, we are able to decouple the indi- 
vidual rate for each transition event from the time series 
(c.f. footnote 42). The results are shown in Fig. 



IV. DISCUSSION 

The diffusion constant of a monomer in our simulations 
is measured to be 1.1 x 10~^ nm^/ps (plots not shown). 
Since the combined binding area for the monomer and 
fibril's ends is about 5 nm^, we expect that the collision 
frequency in our system (with [m] = [/] =7.7 mM) to be 
about 5 X 10~^ ps~^. This is comparable with the rate 
of transition from state 5*0 to state Si^2 observed in our 
simulations (c.f. Fig. [7j), we therefore conclude that the 
initial binding event is well described by a diffusion con- 
trolled reaction. We also note that the initial dissociation 
rate (from state 5'i^2 to state 5*0) is determined to be 
4.3 X 10~^ ps~^, which is comparable to our estimate for 
a- in Eq. 

Besides the transition rate between and 81^2, we 
can see that many of the forward transition rates are of 
the same order of magnitude as the first binding rate. 



FIG. 5: A segment of the data from the simulations, (a) The 
original time series consists of the times when the number 
of D-bond is changed. This is then transformed into a time 
series on the transitions of the set of states {S} (b). The 
procedure of transformation is described in the text. 
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(b) 
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FIG. 6: The numbers of transitions within the set of states 
{S}. Note the ij-entiy denotes the number of transitions from 
state Si-i to Si-i^j. 
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This demonstrates that the elongation process is not 
first-order, but rather dominated by frustrations for the 
monomer to find the correct configuration to become 
fully part of the fibril, i.e., state (c.f. Fig. [1]). We 
note that this qualitative picture stays the same when 
the hydrophobic strength (hydrogen bond strength) are 
varied around the presented values by twenty (ten) per- 
cents above and below. This is the main result of this 
work. 

To have a conceptual feeling for how a rugged energy 
landscape picture would affect the elongation process, we 
look at the work by Zwanzig [sO], which demonstrates 
that: if a particle is diffusing over a ID rugged landscape 
such that the fiuctuation in potential energy is Gaus- 
sian distributed with zero mean and standard deviation 
e, then the motion of the particle can be effectively de- 
scribed by ordinary diffusion with a re-defined diffusion 
constant, I)*, of the form: 



D'' =DeM-{e/kBT)^] 



(7) 



6 



FIG. 7: The transition rates for the set of states {S} con- 
structed from the time series as described in the text. The 
units are in ps~^. 
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FIG. 8: A schematic diagram depicting the scenario where 
the elongation process is viewed as a diffusion process over a 
rugged energy landscape in ID. The linear dimension denotes 
the 'reaction coordinate' and the position xo indicates the 
location when the monomer is first bound to the fibril's end. 
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where D is the original diffusion coefficient. 

Let us now consider the elongation process as a drift- 
diffusion process on a rugged energy landscape (c.f. Fig. 
[8]). Adopting the idea of Zwanzig mentioned above [30], 
we account the ruggedness by redefining the diffusion 
constant as in Eq. ([7j). In other words, the probability 
distribution, of the state of the system (repre- 

sented by the reaction coordinate x) follows the differen- 
tial equation below: 

dtp{x, t) = Ddlp{x, t) - vdxP{x, t) . (8) 

where D is the renormalized diffusion constant that takes 
the ruggedness into account, and v is the drift produced 
by the free energy descent that drives the monomer to 
become fibrilized. The differential equation is supple- 
mented by the boundary condition p(0, t) = p{L^t) = 
where the left boundary depicts monomer detachment 
from the fibril's end and the right boundary depicts com- 
pletion of the fibrilization process (c.f. Fig. [8|). We will 
now assume that the initial condition is a delta function 
located at aL, i.e., p(x, t = 0) = 5{x — aL)^ such that a is 
small. If V is negligible, i.e., when the free energy drive for 
fibrilization is negligible, the ratio of monomers exiting 
at the left boundary (becoming detached) and exiting at 
the right boundary (becoming fibrilized) is proportional 
to a [31] . Using the number of hydrogen bonds again as 
a very crude estimate for the reaction coordinate. For 
the case Abeta peptides, the total number of hydrogen 
bonds is likely to be in the order of 20 [32] so if we take 
the initial location as having one hydrogen bond formed 
with the fibril, a ~ 1/20. In other words, according to 
this diffusion-on-rugged-landscape model, only about 20 
monomers would be recycled before one of them is fibril- 
ized, as compared to the lO'' monomers predicted by the 
two-step model depicted in Eq. ([1]). 

Our models also provides the following experimentally 
relevant insights on the elongation process: 

1. During the period of conversion, the monomer will 
go through a lot of different conformations, and 
there is not a specific conformation that acts as the 
typical conformation before fibrilization. 

2. Since the conversion step is slow, the interac- 



tions between multiple monomers at the fibrils' 
ends should be important, and this propensity 
for monomers' interactions may also serve to pro- 
mote oligomers formation. One would also ex- 
pect that multi-monomer interactions would induce 
more ruggedness into the landscape picture. 

3. Since elongation rates are determined by the form 
of the energy traps, it is expected that the more 
uniform the amyloid-forming peptide's sequence is, 
the slower the elongation rate. This is because 
primary sequence with many identical side chains 
would promote misalignment binding and as a re- 
sult, enhance the ruggedness of the energy land- 
scape. In other words, the complexity of the pri- 
mary sequence may serve as a factor in elongation 
rate prediction (c.f. [H, [m, .38.] ) . 

V. CONCLUSION 

We have studied the elongation process of amyloid fib- 
ril by performing Langevin simulations on a toy model 
of peptides. By projecting the elongation process onto 
a set of discrete states, a rugged energy landscape pic- 
ture emerged, which indicates that monomer-fibril inter- 
action is prolonged in the course of elongation. A crude 
estimate based on this scheme indicates that in the or- 
ders of tens of monomers would have interacted with the 
end of the fibril before a new monomer is fibrilised. Our 
findings also suggest that the complexity of an amyloid 
forming peptide, as measured for instance by how diverse 
the amino- acid compositions are, may serve as a predic- 
tor of the fibril elongation rate. These conclusions can 
be tested with current experimental techniques. 
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